
#define  zprint if (MPIRank==0) print 
!#define zprint print 

module myprecision 
 implicit none  
    integer, parameter         :: sp = KIND(1.0)
    integer, parameter         :: dl = KIND(1.d0)
    integer, parameter         :: pr = sp
    real(pr), parameter        :: fourPI=4.0_pr*3.1415926535897932384626433_pr    
    real(pr), parameter        :: twoPI=2.0_pr*3.1415926535897932384626433_pr    
    real(pr), parameter        :: sqrt_twoPI=sqrt(twoPI)
    real(pr), parameter        :: myee = 2.71828183_pr
end module myprecision 


module module_mychi2 
 use  myprecision
 use  IO
 use  IniFile
 !use  settings
 !use  utils
 !use  ParamDef
 implicit none  

    logical                     :: do_ini0 = .true.
    character(len=100)          :: finname   = "INPUT"  
    character(len=100)          :: fname_exp = "si.Spe"  
    integer                     :: nlines_exp_ini
    real(pr)                    :: dt_exp        !! ns  道宽
    real(pr)                    :: line3_exp
   !real(pr)                    :: fbl_exp       !! ns 分辨率
    integer                     :: ntau_source
    real(pr),allocatable        :: tau_source(:)
    real(pr),allocatable        :: int_source(:)
    real(pr)                    :: intensity_nonsource

#ifdef __5TAU
    integer,parameter           :: ntau_mat = 5  !! 20181021
#endif
#ifdef __4TAU
    integer,parameter           :: ntau_mat = 4  !! 20181021
#endif
#ifdef __3TAU
    integer,parameter           :: ntau_mat = 3  !! 20181021
#endif
#ifdef __2TAU
    integer,parameter           :: ntau_mat = 2  !! 20181021
#endif

    integer                     :: ntau_all
    real(pr),allocatable        :: value_tau(:)  !(1:ntau) != (/6._pr, 0.385_pr, 0.220_pr /)
    real(pr),allocatable        :: intensity(:)  !(1:ntau) != (/0.005_pr , 0.10_pr, 0._pr /)


    integer,allocatable         :: ltpu_exp_ini(:) 
    integer                     :: nleft_chi2 = 6  
    integer                     :: nright_chi2         
    integer                     :: nright_chi2_ini = 400 != 400   !initial


    real(pr),allocatable        :: ltpu_exp(:)          !(1:nleft_chi2+nright_chi2+1) 
   !integer,allocatable         :: ltpu_exp(:)          !(1:nleft_chi2+nright_chi2+1) 
    real(pr),allocatable        :: ltpu_theo(:)         !(1:nleft_chi2+nright_chi2+1) 
    real(pr),allocatable        :: ltpu_d_exp_theo(:)   !(1:nleft_chi2+nright_chi2+1) 

    integer,parameter           :: nbin_perexp=  64     !! 实验数据每道被分割为更多,为了更精细的理论谱  i !! 20181022
    real(pr)                    :: ltpu_background
    integer                     :: ndelta_background = 300
    integer                     :: nstart_background 
    integer                     :: nend_background  
   !integer                     :: nmax_ltpu           !  计算chi2的最高道数
   !integer                     :: nmin_ltpu           !  计算chi2的最低道数

    real(pr),allocatable        :: gauss_array(:)
    integer                     :: gauss_nleft
    integer                     :: gauss_nright
    integer                     :: nleft_theo_t0   ! additional left  range for convolution 
    integer                     :: nright_theo_tf  ! additional right range for convolution 
 

    TYPE                        :: positron_spectrum
        integer                 :: nmin, nmax
        real(pr)                :: xmin, xmax, xdelta
        real(pr)                :: resolution  !! unit : 
        real(pr),allocatable    :: array(:)
        real(pr),allocatable    :: array_c(:)  !! after convolution of array and gauss function
        !integer,allocatable    :: array_c(:)  
    end TYPE 
    TYPE(positron_spectrum)     :: dNofdt !! 光子计数（T） !! /ps
    


 contains 


    subroutine countlines_v2(fname,nlines)
       character(len=100),intent(in)     :: fname
       integer,intent(out)                 :: nlines
       character(len=100)                :: temp
       nlines = 0
       open(1,file=fname) 
       do while(1==1) 
           read(1,*,err=11,end=12) temp
           nlines = nlines + 1
   !        write(*,*) temp
       enddo
11     write(*,*) " error occurs when calling countlines";    stop
12     close(1);return !;write(*,*) "nlines=", nlines 
    end subroutine


    subroutine read_expdata
        character(len=10)       :: temp_10c
        integer                 :: i, temp_int

        integer                 :: myrank, mysize, myerr 

        zprint *, "|read_expdata ..." 
        open(111,file=trim(adjustl( fname_exp ) ) )
            read ( 111, * )  temp_10c
            read ( 111, * )  dt_exp
            read ( 111, * )  temp_10c  !line3_exp
            read ( 111, * )  temp_10c  !fbl_exp
            zprint *, "|dt_exp = ", dt_exp 

            read ( 111, * )  temp_int  !nleft_chi2
            if ( temp_int /= 0 ) then
                nleft_chi2 = temp_int   
                !print*, " set nleft_chi2 = ", nleft_chi2                 
            endif

            read ( 111, * )  ntau_source
            zprint *, "|ntau_source = ", ntau_source 
            allocate(tau_source(1:ntau_source ) )
            allocate(int_source(1:ntau_source ) )

            intensity_nonsource = 1.0
            do i = 1, ntau_source
                read ( 111, FMT=*,err=112,end=113)  tau_source(i), int_source(i)
                zprint *, "|source(", i, ") : tau = ", tau_source(i), " ns, intensity = ", int_source(i)
                intensity_nonsource = intensity_nonsource - int_source(i)
            enddo 
            zprint *, "|intensity_nonsource = ", intensity_nonsource  

            allocate(ltpu_exp_ini(1:nlines_exp_ini-4-1-1-ntau_source ) )
            do i = 1, size(ltpu_exp_ini)
                read ( 111, FMT=*,err=112,end=113)  ltpu_exp_ini(i)
            enddo 
        close(111)
        return 
112  print*, "   error occurs No.111461663 " ;stop
113  print*, "   end No.49813613"; close(123) ; return !;write(*,*) "nlines=", nlines 

    end subroutine 


    subroutine data_preprocessor
        integer                 :: i 
        integer                 :: i_maxloc_exp_ini
        character(len=100)      :: temp_100c

        !character(len=100)      :: fname
        !Type(TIniFile),external :: Ini 
        !character(LEN=*),external  :: InputFile
        !finname = Ini_Read_String('spectrum_filename')
        !print *, Ini_Read_String('spectrum_filename') 
        !fname = Ini_Read_String_File(Ini,'continue_from')
        !print *, "|::", InputFile
        !finname = trim(adjustl( finname ) ) 

        zprint *, "|readin spectrum filename from: "//finname
        open(10,file=trim(adjustl( finname ) ) )
            read ( 10, * )  fname_exp
        close(10)
        zprint *, "|file name of data (exp) : ", trim(adjustl(fname_exp))


        call countlines_v2(fname_exp, nlines_exp_ini)
        zprint *, "|nlines_exp_data = ",  nlines_exp_ini 
        zprint *, "|nbin_perexp       = ", nbin_perexp

        call read_expdata
        nstart_background = int(15.0_pr/dt_exp)  
        nend_background = nstart_background + ndelta_background
        zprint *, "|ndelta_background = ", ndelta_background
        zprint *, "|nstart_background = ", nstart_background
        zprint *, "|nend_background   = ", nend_background

        ntau_all = ntau_source + ntau_mat 
        allocate(value_tau(1:ntau_all ) )
        allocate(intensity(1:ntau_all ) )
        do i = 1, ntau_all
            if ( i .LE. ntau_source ) then
                value_tau(i) = tau_source(i)
                intensity(i) = int_source(i)
            else
                value_tau(i) = 0.0
                intensity(i) = 0.0
            endif
        enddo

        i_maxloc_exp_ini = maxloc(ltpu_exp_ini,1 )
        zprint *, "|max exp value = ", ltpu_exp_ini(i_maxloc_exp_ini), " @ ", i_maxloc_exp_ini

        ltpu_background = sum(ltpu_exp_ini( i_maxloc_exp_ini+nstart_background : & 
                                        i_maxloc_exp_ini+nend_background) ) /(nend_background-nstart_background)
        zprint *, "|ltpu_background = ", ltpu_background                                                    !计算background

        !nmin_ltpu = 1
        nright_chi2=  1
        do while ( ltpu_exp_ini( i_maxloc_exp_ini + nright_chi2) > ltpu_background )
            nright_chi2 = nright_chi2 + 1 
            !print*, "nmax_ltpu =", nmax_ltpu 
        enddo 
        zprint *, "|from max counts, find background counts after ", nright_chi2," channels" 
        !nright_chi2 = floor(nright_chi2 * 1.21_pr / 100.0_pr )*100 +100
        !nright_chi2 = nright_chi2_ini
        nright_chi2 = int(8.0_pr/dt_exp)  
        !print*, " nmin_ltpu    = ", nmin_ltpu 
        !print*, " nmax_ltpu    = ", nmax_ltpu 
        !nright_chi2 = nmax_ltpu !+ 30

        if ( .not. allocated(ltpu_exp) )        allocate( ltpu_exp(1:nleft_chi2+nright_chi2+1) )
        if ( .not. allocated(ltpu_theo) )       allocate( ltpu_theo(1:nleft_chi2+nright_chi2+1) )
        if ( .not. allocated(ltpu_d_exp_theo) ) allocate( ltpu_d_exp_theo(1:nleft_chi2+nright_chi2+1) )
        zprint *, "|set nleft_chi2  = ", nleft_chi2
        zprint *, "|set nright_chi2 = ", nright_chi2
        zprint *, "|total_number_of_chi2: ", nright_chi2+nleft_chi2+1


        !do i = 1, nleft_chi2+nright_chi2+1
        !    ltpu_exp(i) = ltpu_exp_ini(i_maxloc_exp_ini-nleft_chi2-1+i) 
        !enddo
        ltpu_exp = ltpu_exp_ini(i_maxloc_exp_ini-nleft_chi2: i_maxloc_exp_ini+nright_chi2)
        !print*, " ltpu_exp(1:3)", ltpu_exp(1:3)
        !print*, size(ltpu_exp_ini(i_maxloc_exp_ini-nleft_chi2:i_maxloc_exp_ini+nright_chi2))


        dNofdt%xdelta = dt_exp/nbin_perexp 

    end subroutine 
    


    subroutine get_ltpu_theo_ini(temp_ltpu_theo)
        real(pr),intent(out)  ::  temp_ltpu_theo(dNofdt%nmin : dNofdt%nmax )  
        integer                 ::  i , i_ntau, i_start
        real(dl)                ::  dtime1,dtime2

        i_start = (nleft_theo_t0*nbin_perexp + 1)
        !print *, "i_start", i_start

        do i_ntau = 1, ntau_all
        dtime1 = dNofdt%xdelta / value_tau(i_ntau)
        !print *, "i_ntau = ", i_ntau , "dtime1", dtime1

        do i = i_start, dNofdt%nmax                            !!! 2018 nleft

            !if ( mod(i-i_start,1000)<3 ) then 
            !    write(*,'(2I8," ",3G13.6)'),  i_ntau, i, ( i-i_start ) *dNofdt%xdelta /value_tau(i_ntau), dtime2, exp(dtime2)
            !endif
            dtime2 = (i-i_start) * dtime1
            temp_ltpu_theo(i) = temp_ltpu_theo(i) +  intensity(i_ntau)/value_tau(i_ntau) &
                              / exp(dtime2) 
                             !/ exp( ( i-i_start )  * dNofdt%xdelta / value_tau(i_ntau) )
        enddo
        enddo

        !open(661,file="theoini.out.dat")
        ! do i = 1, ( (nleft_chi2+nright_chi2+1+nleft_theo_t0)*nbin_perexp) 
        !  do i = dNofdt%nmax/2.3, dNofdt%nmax/2.3+100
             !write(661,'(f13.8)'), temp_ltpu_theo(i) 
        !     write(*,'(I,", ",f13.8)'), i, temp_ltpu_theo(i) 
        !  end do
        !close(661)
    end subroutine



    function GaussSimpson(AS,BS,NS,FWHM) 
       ! Integral of Gauss_Func form AS to BS  
       IMPLICIT NONE
       real(pr), EXTERNAL           :: FX
       real(pr)                     :: FWHM 
       real(pr)                     :: GaussSimpson,AS,BS
       real(pr)                     :: HS
       real(pr)                     :: zS1,zS2
       real(pr)                     :: SS,fS1,fS2
       INTEGER                     :: NS,I
   
       HS=(BS-AS)/NS
       zS1=AS
       zS2=zS1+HS
       SS=0.0_pr
       fS1=Gauss_Func(zS1,FWHM)
       fS2=Gauss_Func(zS2,FWHM)
       !N=40
       SS = HS*(fS1+4.0_pr*(Gauss_Func( (zS1+zS2)*0.5_pr,FWHM) ) + fS2)/6.0_pr
       DO I = 2, NS
           zS1=zS2
           fS1=fS2
           zS2=zS2+HS
           fS2=Gauss_Func(zS2,FWHM)
           SS = SS + HS*(fS1+4.0_pr*Gauss_Func( (zS1+zS2)*0.5_pr,FWHM) + fS2)/6.0_pr
       END DO
       GaussSimpson = SS
       RETURN
    end function



    function gauss_func(x,fwhm) !!! one dimention 
       IMPLICIT NONE
       real(pr)        :: Gauss_Func,x
       real(pr)        :: FWHM !! 半高全宽 
       real(pr)        :: sigma
       sigma  = FWHM/2.354820047_pr != FWHM/(2._pr*1.1774100236_pr)
       !Gauss_Func = exp(-0.5_pr*x**2.0_pr/sigma**2.0_pr)/sqrt(fourPI/2.0_pr*sigma) 
       !Gauss_Func = exp(-0.5_pr*x**2.0_pr/sigma**2.0_pr) / (sqrt(twoPI)*sigma )
       Gauss_Func = exp(-0.5_pr*(x/sigma)**2.0_pr) !/ ( sqrt_twoPI *sigma )

    end function
    
    

    subroutine get_gauss_array(myFWHM,myxdelta,mynleft,mynright,gauss_array)
        real(pr),intent(in)       :: myFWHM
        real(pr),intent(in)       :: myxdelta
        integer,intent(in)        :: mynleft, mynright
        !real(pr),allocatable      :: get_gauss_array(:) 
        real(pr),intent(out)      :: gauss_array(-mynleft:mynright) 
        real(pr)                  :: myAS, myBS
        integer                   :: i 

        !allocate(get_gauss_array(-mynleft:mynright) ) 
        do i = -mynleft , mynright   
            myAS = i*myxdelta - myxdelta/2._pr 
            myBS = i*myxdelta + myxdelta/2._pr  
            gauss_array(i) = GaussSimpson(myAS,myBS,8,myFWHM)  !! GaussSimpson(myAS,myBS,16,myFWHM)
        enddo

        !gauss_array = gauss_array/sum(gauss_array)
        !if ( sum(get_gauss_array(:)) < 0.999999_pr ) then 
        !    print*, "sum(gauss_array)  = ", sum(get_gauss_array(:)), "< 0.999999_pr, stop program. "
        !    !stop
        !end if 

    end  



#ifdef __5TAU
    function mychi2(p1_nshift, p2_fbl_exp, p3_tau1, p4_intsty1, p5_tau2, p6_intsty2, p7_tau3, p8_intsty3, p9_tau4, p10_intsty4, p11_tau5 ) 
#endif
#ifdef __4TAU
    function mychi2(p1_nshift, p2_fbl_exp, p3_tau1, p4_intsty1, p5_tau2, p6_intsty2, p7_tau3, p8_intsty3, p9_tau4 ) 
#endif
#ifdef __3TAU
    function mychi2(p1_nshift, p2_fbl_exp, p3_tau1, p4_intsty1, p5_tau2, p6_intsty2, p7_tau3)      
#endif
#ifdef __2TAU
    function mychi2(p1_nshift, p2_fbl_exp, p3_tau1, p4_intsty1, p5_tau2)   !! 20201118    
#endif
       !real(pr)                    :: p1_dt_exp
       !real(pr)                    :: p1_dtshift
        real                        :: p1_nshift
        real                        :: p2_fbl_exp
        real                        :: p3_tau1
        real                        :: p4_intsty1
        real                        :: p5_tau2
        real                        :: p6_intsty2
        real                        :: p7_tau3
        real                        :: p8_intsty3
        real                        :: p9_tau4
        real                        :: p10_intsty4
        real                        :: p11_tau5


        !
        real(pr)                    :: mychi2
        real(pr)                    :: mychi2_best
        integer                     :: i_ntau
        integer                     :: i
        integer                     :: id  !! for dNofdt , WofLT
        integer                     :: imaxloc_exp      !实验谱最大数指标 ltpu_exp
        integer                     :: imaxloc_theo_c   !理论谱粗化前的最大数指标 dNofdt%array_c
        integer                     :: iloc_theo_c_max  !理论谱粗化前的最大数指标 dNofdt%array_c
        integer                     :: iloc_theo_c_min  !理论谱粗化前的最小数指标 dNofdt%array_c
        integer                     :: nshift_theo_ini  !粗细理论谱之间变换时的偏移量
        integer                     :: nshift_best
        integer,allocatable,save    :: icountchi2

        real(pr),external           :: ddot
        real(pr),external           :: sdot

        real(pr)                    :: t1, t2

        real, parameter             :: tau_source_1 = 0.385
        !real, parameter             :: intsty_source_1 = 0.1294
        real, parameter             :: intsty_source_1 = 0
        real, parameter             :: tau_source_2 = 2.3134
        !real, parameter             :: intsty_source_2 = 0.0027
        real, parameter             :: intsty_source_2 = 0
        real                        :: intensity_del_source

        !call cpu_time ( t1 )
        !print*, "-------------------------------------------------------------------"
        !print*, "Time0 = ", t1*1000, " ms."

        !print*, "nleft_chi2, nright_chi2 ,nleft_theo_t0 = ", nleft_chi2, nright_chi2, nleft_theo_t0 
        !print*, "nbin_perexp   = ", nbin_perexp
        if( do_ini0 ) then 
            call data_preprocessor 
            do_ini0 = .false.
        endif
        
        !open(112,file="out_ltpu_exp_mod.txt")
        !    write(112,'(1E15.7)') (ltpu_exp(i),i= 1, (nleft_chi2+nright_chi2+1)  ) 
        !close(112)
        !open(112,file="out_ltpu_exp_mod_-bground.txt")
        !    write(112,'(1E15.7)') (ltpu_exp(i) - ltpu_background ,i= 1, (nleft_chi2+nright_chi2+1)  ) 
        !close(112)


        dNofdt%resolution        = p2_fbl_exp
#ifdef __5TAU
        value_tau(ntau_source+1) = p3_tau1
        intensity(ntau_source+1) = intensity_nonsource * p4_intsty1 
        value_tau(ntau_source+2) = p5_tau2
        intensity(ntau_source+2) = intensity_nonsource * p6_intsty2
        value_tau(ntau_source+3) = p7_tau3
        intensity(ntau_source+3) = intensity_nonsource * p8_intsty3
        value_tau(ntau_source+4) = p9_tau4
        intensity(ntau_source+4) = intensity_nonsource * p10_intsty4
        value_tau(ntau_source+5) = p11_tau5
        intensity(ntau_source+5) = intensity_nonsource * (1-p4_intsty1-p6_intsty2-p8_intsty3-p10_intsty4) 
#endif
#ifdef __4TAU
        value_tau(ntau_source+1) = p3_tau1
        intensity(ntau_source+1) = intensity_nonsource * p4_intsty1 
        value_tau(ntau_source+2) = p5_tau2
        intensity(ntau_source+2) = intensity_nonsource * p6_intsty2
        value_tau(ntau_source+3) = p7_tau3
        intensity(ntau_source+3) = intensity_nonsource * p8_intsty3
        value_tau(ntau_source+4) = p9_tau4
        intensity(ntau_source+4) = intensity_nonsource * (1-p4_intsty1-p6_intsty2-p8_intsty3) 
#endif
#ifdef __3TAU
        value_tau(ntau_source+1) = p3_tau1
        intensity(ntau_source+1) = intensity_nonsource * p4_intsty1 
        value_tau(ntau_source+2) = p5_tau2
        intensity(ntau_source+2) = intensity_nonsource * p6_intsty2
        value_tau(ntau_source+3) = p7_tau3
        intensity(ntau_source+3) = intensity_nonsource * (1-p4_intsty1-p6_intsty2) 
#endif
#ifdef __2TAU
        value_tau(ntau_source+1) = p3_tau1
        intensity(ntau_source+1) = intensity_nonsource * p4_intsty1 
        value_tau(ntau_source+2) = p5_tau2
        intensity(ntau_source+2) = intensity_nonsource * (1-p4_intsty1) 
#endif

        
        if ( .not. allocated(icountchi2) ) then
            allocate(icountchi2)
            icountchi2 = 1
        else
            icountchi2 = icountchi2 + 1
        endif


        !if ( abs(sum(intensity) - 1._pr) > 1.E-6_pr )   stop 
        !write(*,'(A,5G13.5)'), " value_tau  = ", value_tau
        !write(*,'(A,5G13.5)'), " intensity  = ", intensity
        !print*,                " dt_exp     = ", p1_dt_exp    
        !print*,                " resolution = ", p2_fbl_exp   
        if ( intensity(ntau_all) < 0._pr) then
            mychi2 = (abs(intensity(ntau_all))+0.01)*999888777666._pr
            zprint *, "|intensity<0: current chi^2 = ", mychi2, ", which has parameters: "
            zprint *, "|intensity<0: p1_nshift = ", p1_nshift, ", p2_fbl_exp = ", p2_fbl_exp
            zprint trim(numcat('(A,',ntau_all) )//'E16.7)' , " |intensity<0: value_tau = ", value_tau(:) 
            zprint trim(numcat('(A,',ntau_all) )//'E16.7)' , " |intensity<0: intensity = ", intensity(:)
            return
        endif


        !call cpu_time ( t1 )
        !print*, "Time1 = ", t1*1000, " ms."

        gauss_nleft = int(3*dNofdt%resolution/dNofdt%xdelta)
        gauss_nright = gauss_nleft 
        if( allocated(gauss_array) ) then
            deallocate(gauss_array)
            allocate( gauss_array( -gauss_nleft : gauss_nright ) )
        else
            allocate( gauss_array( -gauss_nleft : gauss_nright ) )
        endif
        !print*, " gauss_nleft     = ", gauss_nleft
        !print*, " gauss_nright    = ", gauss_nright 

        nleft_theo_t0  = int(gauss_nleft /nbin_perexp) + 7  !! original nleft_theoini
        nright_theo_tf = int(gauss_nright/nbin_perexp) + 7
        !print*, " nleft_theo_t0   = ", nleft_theo_t0
        !print*, " nright_theo_tf  = ", nright_theo_tf

        dNofdt%nmin = 1   !!!  应该不能改变 
        dNofdt%nmax = (nleft_theo_t0+ nleft_chi2+1+nright_chi2 +nright_theo_tf)*nbin_perexp  !!在理论谱零点左边另取nleft_theo_t0个道
        !dNofdt%nmax = (nright_chi2+1+nleft_theo_t0)*nbin_perexp

        dNofdt%xmin = 0.0_pr   !! ns 
        dNofdt%xmax = (dNofdt%nmax-dNofdt%nmin)*dNofdt%xdelta + dNofdt%xmin    !! ns 
        !
        allocate(dNofdt%array(dNofdt%nmin:dNofdt%nmax) )
        allocate(dNofdt%array_c(dNofdt%nmin:dNofdt%nmax) )
        dNofdt%array = 0.0_pr
        dNofdt%array_c = 0.0_pr

        !print*, " dNofdt%nmin     = ", dNofdt%nmin 
        !print*, " dNofdt%nmax     = ", dNofdt%nmax
        !print*, " dNofdt%xdelta   = ", dNofdt%xdelta
        !print*, " dNofdt%xmin     = ", dNofdt%xmin
        !print*, " dNofdt%xmax     = ", dNofdt%xmax

        call get_ltpu_theo_ini(dNofdt%array)
        !open(112,file="out_ltpu_theo_ini.txt" )
        !    write(*,'(10E15.7)') (dNofdt%array(i),i=1, dNofdt%nmax ) 
        !close(112)

        !call cpu_time ( t1 )
        !print*, "Time2 = ", t1*1000, " ms"


        !print*, "gauss_nleft =", gauss_nleft
        !print*, "gauss_nright = ", gauss_nright 
        call get_gauss_array(dNofdt%resolution,dNofdt%xdelta, gauss_nleft, gauss_nright, gauss_array )
        !gauss_array = get_gauss_array(dNofdt%resolution,dNofdt%xdelta, gauss_nleft, gauss_nright)
        !print*, "  gauss_array( gauss_nright )", gauss_array( gauss_nright ) 
        !open(112,file="out_gauss.txt")
        !    write(112,'(1E15.7)') (gauss_array(i),i= -gauss_nleft, gauss_nright ) 
        !close(112)


        !call cpu_time ( t1 )
        !print*, "Time3 = ", t1*1000, " ms"


      if ( 1==1 ) then  !vectorization
        dNofdt%array_c = 0._pr
        do id = dNofdt%nmin+gauss_nleft, dNofdt%nmax-gauss_nright 
          dNofdt%array_c(id) = sum( dNofdt%array(id-gauss_nleft:id+gauss_nright) * gauss_array(:) )
          !if ( pr == sp ) then 
          !  dNofdt%array_c(id) = sdot(gauss_nleft+gauss_nright+1, dNofdt%array(id-gauss_nleft:id+gauss_nright) , 1, gauss_array(-gauss_nleft : gauss_nright), 1)
          !else
          !  dNofdt%array_c(id) = ddot(gauss_nleft+gauss_nright+1, dNofdt%array(id-gauss_nleft:id+gauss_nright) , 1, gauss_array(-gauss_nleft : gauss_nright), 1)
          !endif
        enddo 

      else
        dNofdt%array_c = 0._pr
        do id = dNofdt%nmin, dNofdt%nmax
            do i = id-gauss_nleft, id+gauss_nright
                if (  ( i .GE. (dNofdt%nmin + nleft_theo_t0*nbin_perexp) ) .and. (i .LE. dNofdt%nmax) ) then 
                    dNofdt%array_c(id) = dNofdt%array_c(id) + dNofdt%array(i) * gauss_array(i-id)
                endif
            enddo
        enddo 
        !print*, dNofdt%array_c(dNofdt%nmin+gauss_nleft-2:dNofdt%nmin+gauss_nleft+3)
      end if


        !call cpu_time ( t2 )
        !print*, "Time5 = ",   t2*1000, " ms"
        !print*, "Time5-3   ================== ", (t2-t1)*1000, " ms"


        !!!! 未完成另一种卷积区间: 猜测 imaxloc_theo_c 最小为(nleft_theo_t0*nbin_perexp+1) 
        !!!! 最大为 (nleft_theo_t0*nbin_perexp+gauss_nleft) 
        !!!iloc_theo_c_min  = (nleft_theo_t0*nbin_perexp+1)
        !!!iloc_theo_c_max  = (nleft_theo_t0*nbin_perexp+gauss_nleft)
        !!!dNofdt%array_c = 0._pr
        !!!nshift_theo_ini = 0-  int(nbin_perexp*0.16) ,    nbin_perexp-1  +int(nbin_perexp*0.16 )
        !!!do id = ( imaxloc_theo_c- int((nleft_chi2+0.5)*nbin_perexp)+nshift_theo_ini +1 ) &
        !!!     ,   imaxloc_theo_c- int((nleft_chi2+0.5)*nbin_perexp)+nshift_theo_ini +nbin_perexp*(nleft_chi2+nright_chi2+1) 
        !!!    do i = id-gauss_nleft, id+gauss_nright
        !!!        if (  ( i .GE. dNofdt%nmin ) .and. (i .LE. dNofdt%nmax) ) then 
        !!!            dNofdt%array_c(id) = dNofdt%array_c(id) + dNofdt%array(i) * gauss_array(i-id)
        !!!        endif
        !!!    enddo
        !!!enddo 
        !
        !print*, i, sum(dNofdt%array_c )
        !open(112,file="out_ltpu_theo_conv.txt")
        !    write(112,'(1E15.7)') (dNofdt%array_c(i),i=1, (nleft_chi2+nright+1+nleft_theo_t0)*nbin_perexp ) 
        !close(112)



        mychi2_best = 0.0_pr
        nshift_best = 0
        !nshift_theo_ini = 0 - int(nbin_perexp*0.50), nbin_perexp - 1 + int(nbin_perexp*0.50)
        !p1_nshift = -20
        !nshift_theo_ini = int((p1_dtshift+p1_dt_exp/2) / dNofdt%xdelta) 
        nshift_theo_ini = p1_nshift  !/64.0*nbin_perexp
 
        imaxloc_theo_c      = maxloc( dNofdt%array_c,1 )
        imaxloc_exp         = maxloc( ltpu_exp,1 )
        ltpu_theo           = 0.0_pr 
        do id = 1, (nleft_chi2+nright_chi2+1) 
            ltpu_theo(id) = & 
            sum(dNofdt%array_c(imaxloc_theo_c - int((nleft_chi2+0.5)*nbin_perexp) + nshift_theo_ini + 1 + nbin_perexp*(id-1) & 
                :imaxloc_theo_c- int((nleft_chi2+0.5)*nbin_perexp) + nshift_theo_ini + nbin_perexp*(id)))
            if ( size(ltpu_theo) .EQ. size(ltpu_exp) ) then
                !print *, " theo vs. exp: ", ltpu_theo(id)+ltpu_background , ltpu_exp(id)
            else
                print *, " err : No.4643468741"
                stop
            end if
        enddo 
        !
        ltpu_theo = ltpu_theo / sum(ltpu_theo) * sum(ltpu_exp - ltpu_background)
        ltpu_d_exp_theo = ltpu_theo + ltpu_background - ltpu_exp 

        !open(662,file="ltpu_theo.out.dat")
        !    do i = 1, ( nleft_chi2+nright_chi2+1)
        !        write(662,'(f13.3)'), ltpu_theo(i) 
        !    enddo
        !close(662)

        !call cpu_time ( t1 )
        !print*, "Time6 = ", t1*1000, " ms"

        mychi2 = sum(ltpu_d_exp_theo(1:nleft_chi2+nright_chi2+1)**2/(ltpu_theo(1:nleft_chi2+nright_chi2+1) + ltpu_background))
        !mychi2=sum(ltpu_d_exp_theo(nmin_ltpu:nmax_ltpu)**2/(ltpu_theo(nmin_ltpu:nmax_ltpu)+ltpu_background) )
        !mychi2=sum(ltpu_d_exp_theo(nmin_ltpu:nmax_ltpu)**2/(ltpu_exp(nmin_ltpu:nmax_ltpu) ))
        !!mychi2 = sum( (ltpu_exp(nmin_ltpu:nmax_ltpu)-ltpu_background) **2 / ltpu_background  )
        !stop       
        !print*,  "mychi2 = ", mychi2 

        !call cpu_time ( t1 )
        !print*, "Time7 = ", t1*1000, " ms"

        if ( mychi2_best > mychi2 .or. mychi2_best == 0.0_pr  ) then 
            mychi2_best = mychi2
            nshift_best = nshift_theo_ini 
            !open(112,file="out_ltpu_theo_final.txt")
            !    write(112,'(1E15.7)') (ltpu_theo(i),i= 1, (nleft_chi2+nright_chi2+1)  ) 
            !close(112)
        endif
        !!!!!
        !write(*,*), 'nshift_best = ', nshift_best
        !enddo
        !!!!!!


        if ( mod(icountchi2,2000)==1 ) then
            !do id = 1, (nleft_chi2+nright_chi2+1) 
            !    do i = ( imaxloc_theo_c- int((nleft_chi2+0.5)*nbin_perexp)+nshift_best +1+nbin_perexp*(id-1) ) &
            !     ,   imaxloc_theo_c- int((nleft_chi2+0.5)*nbin_perexp)+nshift_best +nbin_perexp*(id) 
            !       ltpu_theo(id) = ltpu_theo(id) + dNofdt%array_c(i)  
            !    enddo
            !enddo
            !ltpu_theo = ltpu_theo/sum(ltpu_theo)* sum(ltpu_exp-ltpu_background) 
            
            !do i = 1, (nleft_chi2 + 10)
            !    write(*,*), i, ltpu_theo(i) + ltpu_background, ltpu_exp(i)
            !enddo
            !print*, icountchi2, "mychi2_best = ", mychi2_best, "nshift_best = ", nshift_best             
            call cpu_time ( t1 )
            zprint *, "|", icountchi2, " times chi^2, total time cost: ", t1,  " s."
            zprint *, "|current chi^2 = ", mychi2, ", which has parameters: "
            zprint *, "|p1_nshift = ", p1_nshift, ", p2_fbl_exp = ", p2_fbl_exp
            zprint trim(numcat('(A,',ntau_all) )//'E16.7)' , " |value_tau = ", value_tau(:) 
            zprint trim(numcat('(A,',ntau_all) )//'E16.7)' , " |intensity = ", intensity(:)
            zprint *, "|-------------------------------------------------------------------"
            !
            !open(61,file="testSpe.out.dat")
            !do id = 1, (nleft_chi2+nright_chi2+1) 
            !    write(61,'(I4,f13.3,f13.3)'), id, ltpu_theo(id)+ltpu_background , ltpu_exp(id) 
            !enddo
            !close(61)

        end if
        !
        !print*," nshift =",-(0-  int(nbin_perexp*1.16) ) + ( nbin_perexp-1  +int(nbin_perexp*1.16 ))
        !print*, "nmin_ltpu = ", nmin_ltpu 
        !print*, "nmax_ltpu = ", nmax_ltpu 
        !print*, "size(ltpu_exp(nmin_ltpu:nmax_ltpu) ) ", size(ltpu_exp(nmin_ltpu:nmax_ltpu)) 
        !
        !mychi2 = mychi2_best
        !write(*,'( 2E14.7 , 7E14.7 )') dble(size(ltpu_exp(nmin_ltpu:nmax_ltpu))), mychi2,  & 
        !        p2_fbl_exp  ,p3_tau1  ,p4_intsty1  ,p5_tau2  ,p6_intsty2  ,p7_tau3 ,  intensity(ntau_all) 
        !print*,  "mychi2/dof  = ", mychi2_best / size(ltpu_exp(nmin_ltpu:nmax_ltpu)) !/ (nleft_chi2+nright_chi2+1)


        deallocate(dNofdt%array )
        deallocate(dNofdt%array_c )
        deallocate(gauss_array)
    end function

end module module_mychi2 





module CalcLike
 use CMB_Cls
 use cmbtypes
 use cmbdata
 use mpk
 use Random
 use settings
 use ParamDef
 use snovae
 use WeakLen
 use Lya

 !!! ws ad 1line 
 use  myprecision
 use  module_mychi2
 !!!

 implicit none

 logical :: Use_Age_Tophat_Prior = .true.
 logical :: Use_CMB = .true.
 logical :: Use_BBN = .false.
 logical :: Use_Clusters = .false.
  
 integer :: H0_min = 40, H0_max = 100
 real :: Omk_min = -0.3, Omk_max = 0.3
 real :: Use_min_zre = 0
 integer :: Age_min = 10, Age_max = 20
 real :: Temperature  = 1

contains


!  function mychi2(x1,x2,x3)
!    real :: mychi2,x1,x2,x3
!    mychi2 = (x1-1.)**2 + (x2/10.-2.)**2 + (x3/100-3.)**2 
!  end function mychi2


  function GenericLikelihoodFunction(Params) 
    type(ParamSet)  Params 
    real :: GenericLikelihoodFunction
   
  
   !Used when you want to plug in your own CMB-independent likelihood function:
   !set generic_mcmc=.true. in settings.f90, then write function here returning -Ln(Likelihood)
   !Parameter array is Params%P
    !!!!
   ! GenericLikelihoodFunction = greatLike(Params%P)
   ! GenericLikelihoodFunction = LogZero 
   ! call MpiStop('GenericLikelihoodFunction: need to write this function!')

#ifdef __5TAU 
    GenericLikelihoodFunction = mychi2(Params%P(1), Params%P(2), &
                                Params%P(3), Params%P(4), Params%P(5), &
                                Params%P(6), Params%P(7), &
                                Params%P(8), Params%P(9), &
                                Params%P(10), Params%P(11) )
#endif
#ifdef __4TAU 
    GenericLikelihoodFunction = mychi2(Params%P(1), Params%P(2), &
                                Params%P(3), Params%P(4), Params%P(5), &
                                Params%P(6), Params%P(7), &
                                Params%P(8), Params%P(9) )
#endif
#ifdef __3TAU
    GenericLikelihoodFunction = mychi2(Params%P(1), Params%P(2), &
                                Params%P(3), Params%P(4), Params%P(5), &
                                Params%P(6), Params%P(7)) !!! 20181021
#endif
#ifdef __2TAU
    GenericLikelihoodFunction = mychi2(Params%P(1), Params%P(2), &
                                Params%P(3), Params%P(4), Params%P(5) ) !!! 20201118
#endif

    !print*, Params%P(1), Params%P(2), Params%P(3), Params%P(4), Params%P(5)
  end function GenericLikelihoodFunction

  
  function GetLogPrior(CMB, Info) !Get -Ln(Prior)
    real GetLogPrior
    real Age
    Type (CMBParams) CMB
    Type(ParamSetInfo) Info

    GetLogPrior = logZero
 
    if (.not. generic_mcmc) then
     if (CMB%H0 < H0_min .or. CMB%H0 > H0_max) return
     if (CMB%omk < Omk_min .or. CMB%omk > Omk_max .or. CMB%Omv < 0) return
     if (CMB%zre < Use_min_zre) return

     Age = GetAge(CMB, Info)
      !This sets up parameters in CAMB, so do not delete unless you are careful!
 
     if (Use_Age_Tophat_Prior .and. (Age < Age_min .or. Age > Age_max) .or. Age < 0) return
    
    end if
    GetLogPrior = 0
 
  end function GetLogPrior

  function GetLogLike(Params) !Get -Ln(Likelihood)
    type(ParamSet)  Params 
    Type (CMBParams) CMB
    real GetLogLike
    real dum(1,1)
 
    if (any(Params%P > Scales%PMax) .or. any(Params%P < Scales%PMin)) then
       GetLogLike = logZero
        return
    end if

    if (generic_mcmc) then
        GetLogLike = GenericLikelihoodFunction(Params) 
        if (GetLogLike /= LogZero) GetLogLike = GetLogLike/Temperature

    else

     call ParamsToCMBParams(Params%P,CMB)

     GetLogLike = GetLogLikePost(CMB, Params%Info,dum,.false.)
    end if 
   end function GetLogLike

    
  function GetLogLikePost(CMB, Info, inCls, HasCls) 
    real GetLogLikePost
    Type (CMBParams) CMB
    Type(ParamSetInfo) Info
    real, intent(in):: inCls(:,:)
    logical, intent(in) :: HasCls
    real acl(lmax,num_cls_tot)
    integer error

    if (generic_mcmc) stop 'GetLogLikePost: not supported for generic'

    GetLogLikePost  = GetLogPrior(CMB, Info)
    if ( GetLogLikePost >= logZero) then
       GetLogLikePost = logZero
       
    else 
       GetLogLikePost = GetLogLikePost + sum(CMB%nuisance(1:nuisance_params_used)**2)/2
          !Unit Gaussian prior on all nuisance parameters
       if (Use_BBN) GetLogLikePost = GetLogLikePost + (CMB%ombh2 - 0.022)**2/(2*0.002**2) 
          !I'm using increased error bars here
   
       if (Use_CMB .or. Use_LSS) then
          if (HasCls) then
           acl = inCls
           error =0
          else
           call GetCls(CMB, Info, acl, error)
          end if
         if (error /= 0) then
          GetLogLikePost = logZero 
         else
          if (Use_CMB) GetLogLikePost = &
            CMBLnLike(acl, CMB%norm(norm_freq_ix:norm_freq_ix+num_freq_params-1),CMB%nuisance) + GetLogLikePost
          if (Use_mpk) GetLogLikePost = GetLogLikePost + LSSLnLike(CMB, Info%theory)
          if (Use_WeakLen) GetLogLikePost = GetLogLikePost + WeakLenLnLike(CMB, Info%theory)     
          if (Use_Lya) GetLogLikePost = GetLogLikePost +  LSS_Lyalike(CMB, Info%Theory)
          if ( GetLogLikePost >= logZero) then
            GetLogLikePost = logZero
          end if
         end if
         if (Use_SN .and. GetLogLikePost /= logZero ) then
            if (Info%Theory%SN_loglike /= 0) then
             GetLogLikePost = GetLogLikePost + Info%Theory%SN_loglike
            else
             GetLogLikePost = GetLogLikePost + SN_LnLike(CMB)
            end if
               !Assume computed only every time hard parameters change
         end if
         if (Use_BAO .and. GetLogLikePost /= logZero ) then
            if (Info%Theory%BAO_loglike /= 0) then
             GetLogLikePost = GetLogLikePost + Info%Theory%BAO_loglike
            else
             GetLogLikePost = GetLogLikePost + BAO_LnLike(CMB)
            end if
               !Assume computed only every time hard parameters change
         end if
         if (Use_HST .and. GetLogLikePost /= logZero) then
            if (Info%Theory%HST_loglike /= 0) then
             GetLogLikePost = GetLogLikePost + Info%Theory%HST_loglike
            else
             GetLogLikePost = GetLogLikePost + HST_LnLike(CMB)
            end if
           !!Old: GetLogLikePost = GetLogLikePost + (CMB%H0 - 72)**2/(2*8**2)  !HST 
         end if  
     
         else
           if (Use_SN) GetLogLikePost = GetLogLikePost + SN_LnLike(CMB)
           if (Use_BAO) GetLogLikePost = GetLogLikePost + BAO_LnLike(CMB)
           if (Use_HST) GetLogLikePost = GetLogLikePost + HST_LnLike(CMB)           
         end if
         
     
      if (Use_Clusters .and. GetLogLikePost /= LogZero) then
          GetLogLikePost = GetLogLikePost + &
                 (Info%Theory%Sigma_8-0.9)**2/(2*0.05**2)
          stop 'Write your cluster prior in calclike.f90 first!'
      end if

     if (GetLogLikePost /= LogZero) GetLogLikePost = GetLogLikePost/Temperature
   

    end if

  end function GetLogLikePost

end module CalcLike
